Last updated: 2025-11-16
Checks: 5 2
Knit directory: PIPAC_spatial/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown is untracked by Git. To know which version of the R
Markdown file created these results, you’ll want to first commit it to
the Git repo. If you’re still working on the analysis, you can ignore
this warning. When you’re finished, you can run
wflow_publish to commit the R Markdown file and build the
HTML.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20240917) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.
| absolute | relative |
|---|---|
| /home/hnatri/PIPAC_spatial/ | . |
| /home/hnatri/PIPAC_spatial/code/PIPAC_colors_themes.R | code/PIPAC_colors_themes.R |
| /home/hnatri/PIPAC_spatial/code/plot_functions.R | code/plot_functions.R |
| /home/hnatri/PIPAC_spatial/cell_nonimmune_cluster_marker_annotations.tsv | cell_nonimmune_cluster_marker_annotations.tsv |
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 2221fdb. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: analysis/figure/
Ignored: cell_immune_cluster_marker_annotations.tsv
Ignored: cell_nonimmune_cluster_marker_annotations.tsv
Ignored: celltype_markers.tsv
Ignored: code/.Rhistory
Ignored: code/Proximity_analysis/.DS_Store
Ignored: code/Proximity_analysis/data/
Ignored: code/Proximity_analysis/output/
Ignored: code/RSC_latest_EDM_2025-08-06/
Ignored: code/celltype_proximity_tumor.Rout
Ignored: immune_cluster_marker_annotations_2ndpass.tsv
Ignored: immune_prelim_annot_topmarkers.tsv
Ignored: nonimmune_cluster_marker_annotations_2ndpass.tsv
Ignored: nonimmune_prelim_annot_topmarkers.tsv
Untracked files:
Untracked: analysis/finalize_annotations.Rmd
Untracked: analysis/immune_annotation_cell.Rmd
Untracked: analysis/nonimmune_annotation_cell.Rmd
Untracked: code/.RDataTmp
Untracked: code/.RDataTmp1
Untracked: code/celltype_proximity_09232025.R
Untracked: code/celltype_proximity_tumor.R
Untracked: code/denoist_chunk.R
Untracked: code/plot_informative_cluster_markers.R
Untracked: code/slurm.26717417.err
Untracked: code/slurm.26717417.out
Untracked: enrichment_posttreatment_scores_compiled.tsv
Untracked: enrichment_pretreatment_scores_compiled.tsv
Untracked: enrichment_tumor_scores_compiled.tsv
Untracked: output/enrichment_scores_compiled.tsv
Untracked: output/proximity_compiled.tsv
Untracked: tumor_posttreatment_proximity_compiled.tsv
Untracked: tumor_pretreatment_proximity_compiled.tsv
Untracked: tumor_proximity_compiled.tsv
Unstaged changes:
Deleted: Rplots.pdf
Modified: analysis/Xenium_preprocess_wholecell.Rmd
Modified: analysis/add_metadata.Rmd
Modified: analysis/annotation_cell.Rmd
Modified: analysis/arm3_prepost_plots_DEGs.Rmd
Modified: analysis/compare_cellular_nuclear.Rmd
Modified: analysis/compare_cellular_nuclear_denoist_annot.Rmd
Modified: analysis/index.Rmd
Modified: analysis/pca_variance_decomp.Rmd
Modified: analysis/proximity.Rmd
Modified: analysis/splitting_samples.Rmd
Modified: code/PIPAC_colors_themes.R
Modified: code/anndata_to_seurat.R
Modified: code/denoist.R
Modified: code/pairwise_proximity.R
Modified: code/pairwise_proximity_chunk.R
Modified: code/parse_denoist_res.R
Modified: code/plot_functions.R
Modified: code/run_rscript.sh
Modified: code/seurat_to_anndata.R
Modified: code/update_metadata.R
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
There are no past versions. Publish this analysis with
wflow_publish() to start tracking its development.
suppressPackageStartupMessages({
library(workflowr)
library(arrow)
library(dplyr)
library(Seurat)
library(SeuratObject)
library(SeuratDisk)
library(tidyverse)
library(tibble)
library(ggplot2)
library(ggpubr)
library(ggrepel)
library(googlesheets4)
library(workflowr)})
setwd("/home/hnatri/PIPAC_spatial/")
set.seed(9999)
options(scipen = 99999)
options(ggrepel.max.overlaps = Inf)
source("/home/hnatri/PIPAC_spatial/code/PIPAC_colors_themes.R")
source("/home/hnatri/PIPAC_spatial/code/plot_functions.R")
recluster_seurat <- function(seurat_subset, umap_pcs = 10){
seurat_subset <- SCTransform(seurat_subset,
vars.to.regress = c("nCount_RNA", "nFeature_RNA"),
verbose = F)
seurat_subset <- RunPCA(seurat_subset,
npcs = 50,
assay = "SCT",
#features = features,
approx = F,
verbose = F)
integrated_data <- RunUMAP(seurat_subset,
reduction = "pca",
reduction.name = "umap",
dims = 1:umap_pcs,
return.model = TRUE,
verbose = F)
integrated_data <- FindNeighbors(integrated_data,
reduction = "pca",
dims = 1:umap_pcs,
graph.name = c("SCTnn",
"SCTsnn"),
verbose = F)
integrated_data <- FindClusters(integrated_data,
resolution = c(0.1,0.2,0.3,0.5,0.8,1),
graph.name = "SCTsnn",
verbose = F)
return(integrated_data)
}
# Cell type markers by group
gs4_deauth()
markers_by_group <- gs4_get("https://docs.google.com/spreadsheets/d/1w0wrL-5KwNEWi_VpmriN7YmfwpBPlEYX0MKlwI2ZQu8/edit?usp=sharing")
markers_by_group <- read_sheet(markers_by_group, sheet = "Markers by group")
metadata <- gs4_get("https://docs.google.com/spreadsheets/d/1sXXwOreLxjMSUoPt79c6jmaQpluWkaxA5P5HfDsed3I/edit?usp=sharing")
markers <- read_sheet(metadata, sheet = "Markers")
first_pass_annot <- read_sheet(metadata, sheet = "Cell based annotation, non-immune")
seurat_data <- readRDS("/tgen_labs/banovich/PIPAC/Seurat/cell_nonimmune_clustered_NN30_PC20_Seurat.rds")
DefaultAssay(seurat_data) <- "denoist_RNA"
seurat_data <- NormalizeData(seurat_data)
DimPlot(seurat_data,
group.by = "leiden_1.5",
cols = pipac_cell_nonimmune_cluster_col,
reduction = "umap",
raster = T,
label = T) +
coord_fixed(ratio = 1) +
theme_bw() +
NoLegend()
VlnPlot(seurat_data,
features = c("percent.mt", "nCount_cell_RNA", "nFeature_cell_RNA"),
group.by = "leiden_1.5",
ncol = 1,
pt.size = 0) &
theme_bw() &
NoLegend()
seurat_data$leiden_1.5chr <- as.character(seurat_data$leiden_1.5)
seurat_data$leiden_1.5chr <- factor(seurat_data$leiden_1.5chr,
levels = as.character(sort(unique(seurat_data$leiden_1.5))))
DimPlot(seurat_data,
group.by = "leiden_1.5chr",
split.by = "leiden_1.5chr",
cols = pipac_cell_nonimmune_cluster_col,
reduction = "umap",
raster = T,
#label = T,
ncol = 6) +
coord_fixed(ratio = 1) +
theme_bw() +
NoLegend()
Idents(seurat_data) <- seurat_data$leiden_1.5
cluster_markers <- FindAllMarkers(seurat_data,
return.thresh = 0.01,
logfc.threshold = 0.5,
min.pct = 0.20,
verbose = T)
table(cluster_markers$cluster)
7 2 26 15 16 3 21 29 31 12 35 32 25 8 22 6 10 36 24 13
37 43 23 75 25 32 30 50 21 21 140 122 82 75 28 35 25 54 46 38
4 30 33 0 20 23 18 19 9 28 34 17 1 11 27 14 5
44 25 142 34 26 42 86 55 54 27 23 37 33 38 47 52 27
#hist(cluster_markers$avg_log2FC, main = "", xlab = "avg_log2FC", breaks = 100)
#hist(cluster_markers$p_val, main = "", xlab = "p_val", breaks = 100)
#hist(cluster_markers$p_val_adj, main = "", xlab = "p_val_adj", breaks = 100)
top_cluster_markers <- cluster_markers %>%
arrange(dplyr::desc(avg_log2FC)) %>%
group_by(cluster) %>%
dplyr::slice(1:5)
plot_features <- c("PTPRC",
"CD3D", "CD3E", "CD4", "CD8A", # T cells
"STAT4", "STAT3", "TIGIT", "GZMB",
"SELL", "CD19", # B cells
"CD68", "CD44", "MARCO", "APOE", # Macrophages
"C1QB", "C1QBP",
"MUC5AC", "NOTCH3", "RGS5", "MS4A1", "PGA5", "PLVAP", # Lineage markers
"FN1", "DCN", "LUM", # Fibroblasts
"EGR3", "TP53", "JUN", "KIT", # Tumor
"SOX9", "RNF43", "EPCAM")
seurat_data$leiden_1.5 <- factor(seurat_data$leiden_1.5,
levels = sort(unique(seurat_data$leiden_1.5)))
#Idents(seurat_data) <- seurat_data$leiden_1.5
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = plot_features,
cols = c("gray98", "tomato3")) +
RotatedAxis()
FeaturePlot(seurat_data,
features = plot_features,
cols = c("gray98", "tomato3"),
reduction = "umap",
raster = T,
ncol = 5) &
coord_fixed(ratio = 1) &
theme_bw() &
NoLegend()
create_dotplot_heatmap(seurat_object = seurat_data,
plot_features = unique(top_cluster_markers$gene),
group_var = "leiden_1.5",
group_colors = pipac_cell_nonimmune_cluster_col,
column_title = "",
row_km = 5,
col_km = 5,
row.order = NULL,
col.order = NULL)
create_barplot(seurat_data,
plot_var = "Annotation",
group_var = "leiden_1.5",
plot_levels = sort(unique(seurat_data$Annotation)),
group_levels = sort(unique(seurat_data$leiden_1.5)),
plot_colors = pipac_celltype_col,
var_names = c("Annotation", "leiden_1.5"),
legend_title = "")
create_barplot(seurat_data,
group_var = "Annotation",
plot_var = "leiden_1.5",
group_levels = sort(unique(seurat_data$Annotation)),
plot_levels = sort(unique(seurat_data$leiden_1.5)),
plot_colors = pipac_cell_nonimmune_cluster_col,
var_names = c("leiden_1.5", "Annotation"),
legend_title = "")
output_cluster_markers <- cluster_markers %>%
arrange(dplyr::desc(avg_log2FC)) %>%
group_by(cluster) %>%
dplyr::slice(1:30)
output_cluster_markers <- merge(output_cluster_markers, markers, by.x = "gene", by.y = "Gene")
write.table(output_cluster_markers, "/home/hnatri/PIPAC_spatial/cell_nonimmune_cluster_marker_annotations.tsv",
quote = F, row.names = F, sep = "\t")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Cell-cell crosstalk"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Cell-cell crosstalk")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Macrophage signaling"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Macrophage signaling")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "DNA damage and cell faith"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("DNA damage and cell faith")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Functional markers"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Functional markers")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Mesothelial"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Mesothelial")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Macrophages"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Macrophages")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "CAFs"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("CAFs")
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = grep("COL", rownames(seurat_data), value = T),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Collagen")
s.genes <- cc.genes.updated.2019$s.genes
g2m.genes <- cc.genes.updated.2019$g2m.genes
DotPlot(seurat_data,
group.by = "leiden_1.5",
features = c(s.genes, g2m.genes),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("S/G2M genes")
tumor_clusters <- c(7, 8, 25, 32, 33, 35)
stromal_seurat_data <- subset(seurat_data, subset = leiden_1.5 %in% tumor_clusters, invert = T)
#stromal_seurat_data <- recluster_seurat(stromal_seurat_data, vars_to_regress = NULL)
#saveRDS(stromal_seurat_data, "/scratch/hnatri/PIPAC/stromal_seurat_data.rds")
stromal_seurat_data <- readRDS("/tgen_labs/banovich/PIPAC/Seurat/stromal_clustered_NN20_PC20_Seurat.rds")
DefaultAssay(stromal_seurat_data) <- "denoist_RNA"
stromal_seurat_data <- NormalizeData(stromal_seurat_data)
DimPlot(stromal_seurat_data,
group.by = "leiden_1.0",
cols = pipac_cell_nonimmune_cluster_col,
reduction = "umap",
raster = T,
label = T) +
coord_fixed(ratio = 1) +
theme_bw() +
NoLegend()
Idents(stromal_seurat_data) <- stromal_seurat_data$leiden_1.0
cluster_markers <- FindAllMarkers(stromal_seurat_data,
return.thresh = 0.01,
logfc.threshold = 0.5,
min.pct = 0.20,
verbose = T)
top_cluster_markers <- cluster_markers %>%
arrange(dplyr::desc(avg_log2FC)) %>%
group_by(cluster) %>%
dplyr::slice(1:5)
plot_features <- c("PTPRC",
"CD3D", "CD3E", "CD4", "CD8A", # T cells
"STAT4", "STAT3", "TIGIT", "GZMB",
"SELL", "CD19", # B cells
"CD68", "CD44", "MARCO", "APOE", # Macrophages
"C1QB", "C1QBP", "CD163",
"SPP1", "VCAN",
"MUC5AC", "NOTCH3", "RGS5", "MS4A1", # Lineage markers
"PGA5", "PLVAP", "CCL21", "LYVE",
"SPARCL1", # Endo/peri
"FN1", "DCN", "LUM", # Fibroblasts
"COL1A2", "PDPN", "ACTA2", "LMNA",
"EGR3", "TP53", "JUN", "KIT", # Tumor
"SOX9", "RNF43", "EPCAM", "PECAM", "CEACAM")
stromal_seurat_data$leiden_1.0 <- factor(stromal_seurat_data$leiden_1.0,
levels = sort(unique(stromal_seurat_data$leiden_1.0)))
Idents(stromal_seurat_data) <- stromal_seurat_data$leiden_1.0
DotPlot(stromal_seurat_data,
group.by = "leiden_1.0",
features = plot_features,
cols = c("gray98", "tomato3")) +
RotatedAxis()
FeaturePlot(stromal_seurat_data,
features = plot_features,
cols = c("gray98", "tomato3"),
reduction = "umap",
raster = T,
ncol = 5) &
coord_fixed(ratio = 1) &
theme_bw() &
NoLegend()
create_dotplot_heatmap(seurat_object = stromal_seurat_data,
plot_features = unique(top_cluster_markers$gene),
group_var = "leiden_1.0",
group_colors = pipac_cell_nonimmune_cluster_col,
column_title = "",
row_km = 5,
col_km = 5,
row.order = NULL,
col.order = NULL)
DotPlot(stromal_seurat_data,
group.by = "leiden_1.0",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Cell-cell crosstalk"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Cell-cell crosstalk")
DotPlot(seurat_data,
group.by = "leiden_1.0",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "DNA damage and cell faith"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("DNA damage and cell faith")
DotPlot(stromal_seurat_data,
group.by = "leiden_1.0",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Functional markers"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Functional markers")
DotPlot(stromal_seurat_data,
group.by = "leiden_1.0",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "Mesothelial"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Mesothelial")
DotPlot(stromal_seurat_data,
group.by = "leiden_1.0",
features = unique(markers_by_group[which(markers_by_group$annotation1 == "CAFs"),]$feature),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("CAFs")
DotPlot(stromal_seurat_data,
group.by = "leiden_1.0",
features = grep("COL", rownames(seurat_data), value = T),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("Collagen")
s.genes <- cc.genes.updated.2019$s.genes
g2m.genes <- cc.genes.updated.2019$g2m.genes
DotPlot(stromal_seurat_data,
group.by = "leiden_1.0",
features = c(s.genes, g2m.genes),
cols = c("azure", "tomato3")) +
RotatedAxis() +
ggtitle("S/G2M genes")
pipac_celltype_col <- c(pipac_celltype_col, c("NA" = "lightgray"))
stromal_seurat_data$Annotation <- ifelse(is.na(stromal_seurat_data$Annotation), "NA", stromal_seurat_data$Annotation)
create_barplot(seurat_data,
group_var = "Annotation",
plot_var = "leiden_0.5",
group_levels = sort(unique(seurat_data$Annotation)),
plot_levels = sort(unique(seurat_data$leiden_1.0)),
plot_colors = pipac_cluster_20_col,
var_names = c("leiden_0.5", "Annotation"),
legend_title = "")
create_barplot(seurat_data,
plot_var = "Annotation",
group_var = "leiden_1.0",
plot_levels = sort(unique(seurat_data$Annotation)),
group_levels = sort(unique(seurat_data$leiden_1.0)),
plot_colors = pipac_celltype_col,
var_names = c("Annotation", "leiden_0.5"),
legend_title = "")
create_barplot(stromal_seurat_data,
group_var = "Annotation",
plot_var = "leiden_0.5",
group_levels = sort(unique(stromal_seurat_data$Annotation)),
plot_levels = sort(unique(stromal_seurat_data$leiden_1.0)),
plot_colors = pipac_cluster_20_col,
var_names = c("leiden_0.5", "Annotation"),
legend_title = "")
create_barplot(stromal_seurat_data,
plot_var = "Annotation",
group_var = "leiden_1.0",
plot_levels = sort(unique(stromal_seurat_data$Annotation)),
group_levels = sort(unique(stromal_seurat_data$leiden_1.0)),
plot_colors = pipac_celltype_col,
var_names = c("Annotation", "leiden_0.5"),
legend_title = "")
stromal_recluster_subset <- subset(stromal_seurat_data, subset = leiden_1.0 %in% c(0, 15, 23, 26))
stromal_recluster_subset <- recluster_seurat(stromal_recluster_subset, umap_pcs = 10)
DimPlot(stromal_recluster_subset,
group.by = "SCTsnn_res.0.8",
label = T) +
NoLegend()
stromal_recluster_subset$SCTsnn_res.0.8 <- factor(stromal_recluster_subset$SCTsnn_res.0.8,
levels = sort(unique(stromal_recluster_subset$SCTsnn_res.0.8)))
Idents(stromal_recluster_subset) <- stromal_recluster_subset$SCTsnn_res.0.8
DotPlot(stromal_recluster_subset,
#group.by = "leiden_1.0",
features = plot_features,
cols = c("gray89", "tomato3")) +
RotatedAxis()
Idents(stromal_recluster_subset) <- stromal_recluster_subset$SCTsnn_res.0.8
stromal_cluster_markers <- FindAllMarkers(stromal_recluster_subset,
return.thresh = 0.01,
logfc.threshold = 0.5,
min.pct = 0.20,
verbose = T)
table(stromal_cluster_markers$cluster)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
128 73 37 69 97 34 75 66 121 86 39 73 119 101 73 104 93 71 60 79
20
113
stromal_top_cluster_markers <- stromal_cluster_markers %>%
arrange(dplyr::desc(avg_log2FC)) %>%
group_by(cluster) %>%
dplyr::slice(1:5)
create_dotplot_heatmap(seurat_object = stromal_recluster_subset,
plot_features = unique(stromal_top_cluster_markers$gene),
group_var = "SCTsnn_res.0.8",
group_colors = pipac_cell_nonimmune_cluster_col,
column_title = "",
row_km = 5,
col_km = 5,
row.order = NULL,
col.order = NULL)
FeaturePlot(subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 10),
features = c("DCN", "PLVAP"),
blend = T) &
theme_bw()
cl10 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 10)
cl10 <- recluster_seurat(cl10, umap_pcs = 10)
FeaturePlot(cl10,
features = c("DCN", "PLVAP"),
blend = T) &
theme_bw()
FeaturePlot(cl10,
features = c("DCN", "PLVAP")) &
theme_bw()
ncol(cl10)
[1] 1102
table(cl10$SCTsnn_res.0.8)
0 1 2 3 4 5 6
281 183 183 146 135 106 68
DimPlot(cl10,
group.by = "SCTsnn_res.0.8",
reduction = "umap",
label = T) +
coord_fixed(ratio = 1) +
theme_bw() +
NoLegend() +
ggtitle("Cluster 10")
cl10$SCTsnn_res.0.8 <- factor(cl10$SCTsnn_res.0.8,
levels = sort(unique(cl10$SCTsnn_res.0.8)))
Idents(cl10) <- cl10$SCTsnn_res.0.8
DotPlot(cl10,
group.by = "SCTsnn_res.0.8",
features = plot_features,
cols = c("gray98", "tomato3")) +
RotatedAxis()
Idents(cl10) <- cl10$SCTsnn_res.0.8
cl10_markers <- FindAllMarkers(cl10,
return.thresh = 0.01,
logfc.threshold = 0.5,
min.pct = 0.20,
verbose = F)
table(cl10_markers$cluster)
0 1 2 3 4 5 6
32 37 23 31 43 43 41
cl10_markers %>%
arrange(dplyr::desc(avg_log2FC)) %>%
group_by(cluster) %>%
dplyr::slice(1:5)
# A tibble: 35 × 7
# Groups: cluster [7]
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
1 1.09e-57 2.20 0.808 0.345 4.90e-55 0 MGP
2 6.47e-58 1.94 0.783 0.262 2.90e-55 0 DCN
3 3.78e-55 1.85 0.786 0.253 1.69e-52 0 CXCL14
4 3.12e-18 1.60 0.363 0.128 1.40e-15 0 GPNMB
5 7.72e-22 1.55 0.491 0.205 3.46e-19 0 LUM
6 3.98e-66 3.70 0.541 0.064 1.78e-63 1 CCL14
7 1.33e-29 2.83 0.475 0.144 5.97e-27 1 POSTN
8 6.37e-66 2.53 0.847 0.258 2.85e-63 1 PLVAP
9 4.47e-11 1.96 0.224 0.07 2.00e- 8 1 ENTPD1
10 2.02e-29 1.92 0.557 0.171 9.05e-27 1 FLT1
# ℹ 25 more rows
nos <- c()
fibro1 <- c() # FN1 high
fibro2 <- c(0, 4) # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c(3) # NOTCH3, RGS5
endo <- c(1, 2, 5) # PLVAP
myel <- c(6)
cl10$annotation <- ifelse(cl10$SCTsnn_res.0.8 %in% peri, "Peri_from_stromal",
ifelse(cl10$SCTsnn_res.0.8 %in% endo, "Endo_from_stromal",
ifelse(cl10$SCTsnn_res.0.8 %in% fibro2, "F2_from_stromal",
ifelse(cl10$SCTsnn_res.0.8 %in% myel, "Myel_from_stromal", NA))))
cl10$annotation <- paste0(cl10$SCTsnn_res.0.8, "_", cl10$annotation)
unique(cl10$annotation)
[1] "0_F2_from_stromal" "5_Endo_from_stromal" "1_Endo_from_stromal"
[4] "4_F2_from_stromal" "2_Endo_from_stromal" "6_Myel_from_stromal"
[7] "3_Peri_from_stromal"
FeaturePlot(subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 11),
features = c("FN1", "CEACAM6"),
blend = T) &
theme_bw()
reccl11 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 11)
reccl11 <- recluster_seurat(reccl11, umap_pcs = 10)
FeaturePlot(reccl11,
features = c("DCN", "MARCO"),
blend = T) &
theme_bw()
FeaturePlot(reccl11,
features = c("DCN", "MARCO")) &
theme_bw()
ncol(reccl11)
[1] 1100
table(reccl11$SCTsnn_res.0.8)
0 1 2 3 4 5 6 7 8
240 205 199 97 87 85 75 69 43
DimPlot(reccl11,
group.by = "SCTsnn_res.0.8",
reduction = "umap",
label = T) +
coord_fixed(ratio = 1) +
theme_bw() +
NoLegend() +
ggtitle("Cluster 11")
reccl11$SCTsnn_res.0.8 <- factor(reccl11$SCTsnn_res.0.8,
levels = sort(unique(reccl11$SCTsnn_res.0.8)))
Idents(reccl11) <- reccl11$SCTsnn_res.0.8
DotPlot(reccl11,
group.by = "SCTsnn_res.0.8",
features = plot_features,
cols = c("gray98", "tomato3")) +
RotatedAxis()
nos <- c(3)
fibro1 <- c(0, 1, 2, 4, 7, 8) # FN1 high
fibro2 <- c() # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c() # NOTCH3, RGS5
endo <- c() # PLVAP
myel <- c(5, 6)
reccl11$annotation <- ifelse(reccl11$SCTsnn_res.0.8 %in% nos, "NOS_from_stromal",
ifelse(reccl11$SCTsnn_res.0.8 %in% fibro1, "F1_from_stromal",
ifelse(reccl11$SCTsnn_res.0.8 %in% myel, "Myel_from_stromal", NA)))
reccl11$annotation <- paste0(reccl11$SCTsnn_res.0.8, "_", reccl11$annotation)
unique(reccl11$annotation)
[1] "2_F1_from_stromal" "0_F1_from_stromal" "8_F1_from_stromal"
[4] "3_NOS_from_stromal" "1_F1_from_stromal" "4_F1_from_stromal"
[7] "7_F1_from_stromal" "5_Myel_from_stromal" "6_Myel_from_stromal"
FeaturePlot(subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 13),
features = c("COL1A2", "APOE"),
blend = T) &
theme_bw()
reccl13 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 13)
reccl13 <- recluster_seurat(reccl13, umap_pcs = 10)
FeaturePlot(reccl13,
features = c("COL1A2", "APOE"),
blend = T) &
theme_bw()
FeaturePlot(reccl13,
features = c("COL1A2", "APOE")) &
theme_bw()
ncol(reccl13)
[1] 910
table(reccl13$SCTsnn_res.0.8)
0 1 2 3 4 5 6
216 170 134 104 104 94 88
DimPlot(reccl13,
group.by = "SCTsnn_res.0.8",
reduction = "umap",
label = T) +
coord_fixed(ratio = 1) +
theme_bw() +
NoLegend() +
ggtitle("Cluster 13")
reccl13$SCTsnn_res.0.8 <- factor(reccl13$SCTsnn_res.0.8,
levels = sort(unique(reccl13$SCTsnn_res.0.8)))
Idents(reccl13) <- reccl13$SCTsnn_res.0.8
DotPlot(reccl13,
group.by = "SCTsnn_res.0.8",
features = plot_features,
cols = c("gray98", "tomato3")) +
RotatedAxis()
nos <- c()
fibro1 <- c(1, 2) # FN1 high
fibro2 <- c() # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c() # NOTCH3, RGS5
endo <- c() # PLVAP
myel <- c(0, 3, 4, 5, 6)
reccl13$annotation <- ifelse(reccl13$SCTsnn_res.0.8 %in% fibro1, "F1_from_stromal",
ifelse(reccl13$SCTsnn_res.0.8 %in% myel, "Myel_from_stromal", NA))
reccl13$annotation <- paste0(reccl13$SCTsnn_res.0.8, "_", reccl13$annotation)
unique(reccl13$annotation)
[1] "0_Myel_from_stromal" "2_F1_from_stromal" "5_Myel_from_stromal"
[4] "3_Myel_from_stromal" "4_Myel_from_stromal" "1_F1_from_stromal"
[7] "6_Myel_from_stromal"
reccl20 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 20)
reccl20 <- recluster_seurat(reccl20, umap_pcs = 10)
FeaturePlot(reccl20,
features = c("COL1A2", "PLVAP"),
blend = T) &
theme_bw()
FeaturePlot(reccl20,
features = c("COL1A2", "PLVAP")) &
theme_bw()
ncol(reccl20)
[1] 247
table(reccl20$SCTsnn_res.0.8)
0 1 2 3
80 78 48 41
DimPlot(reccl20,
group.by = "SCTsnn_res.0.8",
reduction = "umap",
label = T) +
coord_fixed(ratio = 1) +
theme_bw() +
NoLegend() +
ggtitle("Cluster 20")
reccl20$SCTsnn_res.0.8 <- factor(reccl20$SCTsnn_res.0.8,
levels = sort(unique(reccl20$SCTsnn_res.0.8)))
Idents(reccl20) <- reccl20$SCTsnn_res.0.8
DotPlot(reccl20,
group.by = "SCTsnn_res.0.8",
features = plot_features,
cols = c("gray98", "tomato3")) +
RotatedAxis()
nos <- c()
fibro1 <- c(0, 1, 3) # FN1 high
fibro2 <- c() # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c(2) # NOTCH3, RGS5
endo <- c() # PLVAP
myel <- c()
reccl20$annotation <- ifelse(reccl20$SCTsnn_res.0.8 %in% fibro1, "F1_from_stromal",
ifelse(reccl20$SCTsnn_res.0.8 %in% peri, "Peri_from_stromal", NA))
reccl20$annotation <- paste0(reccl20$SCTsnn_res.0.8, "_", reccl20$annotation)
# Merging
seurat_list <- ls(pattern="reccl")
seurat_list <- str_sort(seurat_list, numeric = TRUE)
seurat_list <- do.call("list", mget(seurat_list))
#saveRDS(seurat_list, "/scratch/hnatri/PIPAC/seurat_list.Rds")
seurat_list <- readRDS("/scratch/hnatri/PIPAC/seurat_list.Rds")
seurat_list <- lapply(names(seurat_list), function(xx){
#message(xx)
cll <- xx
xx <- seurat_list[[xx]]
xx[["denoist_RNA"]] <- as(object = xx[["denoist_RNA"]], Class = "Assay5")
DefaultAssay(xx) <- "denoist_RNA"
xx[["RNA"]] <- NULL
xx[["nucleus_RNA"]] <- NULL
xx[["SCT"]] <- NULL
xx$annotation <- paste0(cll, "_", xx$annotation)
xx@meta.data <- xx@meta.data[,c("cell_id", "annotation")]
xx
})
seurat_merged <- merge(x = seurat_list[[1]], y = seurat_list[2:3])
saveRDS(seurat_merged, "/scratch/hnatri/PIPAC/seurat_merged.rds")
#seurat_merged <- readRDS("/scratch/hnatri/PIPAC/seurat_merged.rds")
# Adding most granular annotations to the stromal reclustered subset object
stromal_recluster_subset$stromal_subcluster_annot <- mapvalues(x = colnames(stromal_recluster_subset),
from = colnames(seurat_merged),
to = seurat_merged$annotation)
stromal_recluster_subset@meta.data[-which(colnames(stromal_recluster_subset) %in% colnames(seurat_merged)),]$stromal_subcluster_annot <- NA
# Adding annotations to the complete stromal subset
stromal_ambig_annot <- read_sheet(metadata, sheet = "Stromal ambiguous recluster")
#stromal_ambig_annot <- stromal_ambig_annot[-which(is.na(stromal_ambig_annot$annotation)),]
#stromal_ambig_annot$annotation <- paste0(stromal_ambig_annot$SCTsnn_res.0.8, "_", stromal_ambig_annot$annotation)
stromal_seurat_data$stromal_subcluster_annot <- mapvalues(x = colnames(stromal_seurat_data),
from = colnames(stromal_recluster_subset),
to = stromal_recluster_subset$stromal_subcluster_annot)
stromal_seurat_data@meta.data[-which(colnames(stromal_seurat_data) %in% colnames(stromal_recluster_subset)),]$stromal_subcluster_annot <- NA
stromal_recluster_subset$stromal_ambig_annot <- mapvalues(x = stromal_recluster_subset$SCTsnn_res.0.8,
from = stromal_ambig_annot$SCTsnn_res.0.8,
to = as.character(stromal_ambig_annot$annotation))
#stromal_recluster_subset@meta.data[-which(stromal_recluster_subset$SCTsnn_res.0.8 %in% stromal_ambig_annot$SCTsnn_res.0.8),]$stromal_ambig_annot <- NA
# Adding stromal annotations to the full non-immune object
seurat_data$stromal_ambig_annot <- mapvalues(x = colnames(seurat_data),
from = colnames(stromal_recluster_subset),
to = as.character(stromal_recluster_subset$stromal_ambig_annot))
seurat_data@meta.data[-which(colnames(seurat_data) %in% colnames(stromal_recluster_subset)),]$stromal_ambig_annot <- NA
seurat_data$stromal_subcluster_annot <- mapvalues(x = colnames(seurat_data),
from = colnames(stromal_recluster_subset),
to = stromal_recluster_subset$stromal_subcluster_annot)
seurat_data@meta.data[-which(colnames(seurat_data) %in% colnames(seurat_merged)),]$stromal_subcluster_annot <- NA
# Adding labels for immune clusters in the whole stromal subset
stromal_seurat_data$stromal_annot <- ifelse(stromal_seurat_data$leiden_1.0 == 1, "Immune1_from_stromal",
ifelse(stromal_seurat_data$leiden_1.0 == 11, "Immune2_from_stromal",
ifelse(stromal_seurat_data$leiden_1.0 == 25, "Immune3_from_stromal", NA)))
seurat_data$stromal_annot <- mapvalues(x = colnames(seurat_data),
from = colnames(stromal_seurat_data),
to = stromal_seurat_data$stromal_annot)
seurat_data@meta.data[-which(colnames(seurat_data) %in% colnames(stromal_seurat_data)),]$stromal_annot <- NA
# Adding annotation for the whole non-immune subset
# first_pass_annot
seurat_data$nonimmune_annotation <- mapvalues(x = seurat_data$leiden_1.5,
from = first_pass_annot$leiden_1.5,
to = first_pass_annot$annotation)
seurat_data@meta.data[-which(seurat_data$leiden_1.5 %in% stromal_seurat_data$leiden_1.5),]$nonimmune_annotation <- NA
seurat_data$nonimmune_annotation <- as.character(seurat_data$nonimmune_annotation)
# Combining
#saveRDS(seurat_data, "/scratch/hnatri/PIPAC/nonimmune_annotated.rds")
#seurat_data <- readRDS("/scratch/hnatri/PIPAC/nonimmune_annotated.rds")
seurat_data$combined_annotation <- ifelse(is.na(seurat_data$stromal_subcluster_annot), seurat_data$stromal_ambig_annot, seurat_data$stromal_subcluster_annot)
seurat_data$combined_annotation <- ifelse(is.na(seurat_data$combined_annotation), as.character(seurat_data$stromal_annot), as.character(seurat_data$combined_annotation))
seurat_data$combined_annotation <- ifelse(is.na(seurat_data$combined_annotation), seurat_data$nonimmune_annotation, seurat_data$combined_annotation)
#seurat_data$combined_annotation <- ifelse(is.na(seurat_data$combined_annotation), seurat_data$annotation, #NA)
saveRDS(seurat_data, "/scratch/hnatri/PIPAC/nonimmune_annotated.rds")
# Rscript -e "rmarkdown::render('nonimmune_annotation_cell.Rmd')"
# mv *html /home/hnatri/PIPAC_spatial/docs/
sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] ComplexHeatmap_2.18.0 viridis_0.6.3 viridisLite_0.4.2
[4] circlize_0.4.15 plyr_1.8.8 RColorBrewer_1.1-3
[7] googlesheets4_1.1.0 ggrepel_0.9.3 ggpubr_0.6.0
[10] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
[13] purrr_1.0.1 readr_2.1.4 tidyr_1.3.0
[16] tibble_3.2.1 ggplot2_3.4.2 tidyverse_2.0.0
[19] SeuratDisk_0.0.0.9021 Seurat_5.0.1 SeuratObject_5.0.1
[22] sp_1.6-1 dplyr_1.1.2 arrow_21.0.0.1
[25] workflowr_1.7.1
loaded via a namespace (and not attached):
[1] fs_1.6.2 matrixStats_1.0.0
[3] spatstat.sparse_3.0-1 bitops_1.0-7
[5] httr_1.4.6 doParallel_1.0.17
[7] tools_4.3.0 sctransform_0.4.1
[9] backports_1.4.1 utf8_1.2.3
[11] R6_2.5.1 lazyeval_0.2.2
[13] uwot_0.1.14 GetoptLong_1.0.5
[15] withr_2.5.0 gridExtra_2.3
[17] progressr_0.13.0 cli_3.6.1
[19] Biobase_2.62.0 Cairo_1.6-0
[21] spatstat.explore_3.2-1 fastDummies_1.7.3
[23] labeling_0.4.2 sass_0.4.6
[25] spatstat.data_3.0-1 ggridges_0.5.4
[27] pbapply_1.7-0 parallelly_1.36.0
[29] limma_3.58.1 rstudioapi_0.14
[31] generics_0.1.3 shape_1.4.6
[33] ica_1.0-3 spatstat.random_3.1-5
[35] car_3.1-2 Matrix_1.6-5
[37] ggbeeswarm_0.7.2 fansi_1.0.4
[39] S4Vectors_0.40.2 abind_1.4-5
[41] lifecycle_1.0.3 whisker_0.4.1
[43] yaml_2.3.7 carData_3.0-5
[45] SummarizedExperiment_1.32.0 SparseArray_1.2.3
[47] Rtsne_0.16 glmGamPoi_1.14.3
[49] promises_1.2.0.1 crayon_1.5.2
[51] miniUI_0.1.1.1 lattice_0.21-8
[53] cowplot_1.1.1 magick_2.7.4
[55] pillar_1.9.0 knitr_1.43
[57] GenomicRanges_1.54.1 rjson_0.2.21
[59] future.apply_1.11.0 codetools_0.2-19
[61] leiden_0.4.3 glue_1.6.2
[63] getPass_0.2-4 data.table_1.14.8
[65] vctrs_0.6.2 png_0.1-8
[67] spam_2.9-1 cellranger_1.1.0
[69] gtable_0.3.3 assertthat_0.2.1
[71] cachem_1.0.8 xfun_0.39
[73] S4Arrays_1.2.0 mime_0.12
[75] survival_3.5-5 gargle_1.4.0
[77] iterators_1.0.14 statmod_1.5.0
[79] ellipsis_0.3.2 fitdistrplus_1.1-11
[81] ROCR_1.0-11 nlme_3.1-162
[83] bit64_4.0.5 RcppAnnoy_0.0.20
[85] GenomeInfoDb_1.38.5 rprojroot_2.0.3
[87] bslib_0.4.2 irlba_2.3.5.1
[89] vipor_0.4.5 KernSmooth_2.23-21
[91] colorspace_2.1-0 BiocGenerics_0.48.1
[93] ggrastr_1.0.2 tidyselect_1.2.0
[95] processx_3.8.1 bit_4.0.5
[97] compiler_4.3.0 curl_5.0.0
[99] git2r_0.32.0 hdf5r_1.3.8
[101] DelayedArray_0.28.0 plotly_4.10.2
[103] scales_1.2.1 lmtest_0.9-40
[105] callr_3.7.3 digest_0.6.31
[107] goftest_1.2-3 spatstat.utils_3.0-3
[109] presto_1.0.0 rmarkdown_2.22
[111] XVector_0.42.0 htmltools_0.5.5
[113] pkgconfig_2.0.3 sparseMatrixStats_1.14.0
[115] MatrixGenerics_1.14.0 highr_0.10
[117] fastmap_1.1.1 rlang_1.1.1
[119] GlobalOptions_0.1.2 htmlwidgets_1.6.2
[121] DelayedMatrixStats_1.24.0 shiny_1.7.4
[123] farver_2.1.1 jquerylib_0.1.4
[125] zoo_1.8-12 jsonlite_1.8.5
[127] RCurl_1.98-1.12 magrittr_2.0.3
[129] GenomeInfoDbData_1.2.11 dotCall64_1.0-2
[131] patchwork_1.1.2 munsell_0.5.0
[133] Rcpp_1.0.10 reticulate_1.29
[135] stringi_1.7.12 zlibbioc_1.48.0
[137] MASS_7.3-60 parallel_4.3.0
[139] listenv_0.9.0 deldir_1.0-9
[141] splines_4.3.0 tensor_1.5
[143] hms_1.1.3 ps_1.7.5
[145] igraph_1.4.3 spatstat.geom_3.2-1
[147] ggsignif_0.6.4 RcppHNSW_0.5.0
[149] reshape2_1.4.4 stats4_4.3.0
[151] evaluate_0.21 tzdb_0.4.0
[153] foreach_1.5.2 httpuv_1.6.11
[155] RANN_2.6.1 polyclip_1.10-4
[157] future_1.32.0 clue_0.3-64
[159] scattermore_1.2 broom_1.0.4
[161] xtable_1.8-4 RSpectra_0.16-1
[163] rstatix_0.7.2 later_1.3.1
[165] googledrive_2.1.0 beeswarm_0.4.0
[167] IRanges_2.36.0 cluster_2.1.4
[169] timechange_0.2.0 globals_0.16.2